using Plots
using Images
using ImageDraw
using VideoIO
using Base64
using CSV
using DataFrames
using DifferentialEquations.OrdinaryDiffEq
using Optim
using StatsBase
include("scripts/DisplayMovie.jl")
include("scripts/ImageTracking.jl")
Main.ImageTracking
filename = "img/pendulo_70cm_1_reduzido_curto.mp4"
vd = VideoIO.load(filename)
typeof(copy(vd[1]))
Matrix{RGB{N0f8}} (alias for Array{RGB{Normed{UInt8, 8}}, 2})
println("Duração do vÃdeo: $(VideoIO.get_duration(filename)) s")
println("Número de quadros: $(length(vd))")
println("Número médio de quadros por segundo: $(length(vd)/VideoIO.get_duration(filename))")
Duração do vÃdeo: 5.086 s Número de quadros: 150 Número médio de quadros por segundo: 29.49272512780181
# video completo
filename_completo = "img/pendulo_70cm_1_reduzido.mov"
vd_completo = VideoIO.load(filename_completo)
println("""Duração do vÃdeo: $(VideoIO.get_duration(filename_completo)) s""")
println("Número de quadros: $(length(vd_completo))")
println("""Número médio de quadros por segundo: $(length(vd_completo)/VideoIO.get_duration(filename_completo))""")
Duração do vÃdeo: 16.45 s Número de quadros: 493 Número médio de quadros por segundo: 29.969604863221885
vd[1]
DisplayMovie.display_movie("img/pendulo_70cm_1_reduzido_curto.mp4")
background = convert.(RGB{Float16}, vd[1])
convert.(RGB{Float16}, vd[54])
convert.(Gray{Float16}, abs.(convert.(RGB{Float16}, vd[54]) - background))
threshold = 0.1
mask = convert.(
Gray{Float16},
abs.(
convert.(RGB{Float16}, vd[54]) -
background
)
) .> threshold
480×854 BitMatrix: 0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 … 1 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 … 1 1 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 1 1 1 1 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
convert.(Gray{Float16}, mask)
labeled_components = label_components(mask)
480×854 Matrix{Int64}:
0 0 0 0 0 0 0 0 0 0 0 0 0 … 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 … 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 … 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249
⋮ ⋮ ⋮ ⋱ ⋮
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
@which label_components(mask)
?label_components
search: label_components label_components! labeled_components
label = label_components(tf, [connectivity])
label = label_components(tf, [region])
Find the connected components in a binary array tf. There are two forms that connectivity can take:
tf, of size 1 or 3 along each dimension. Each entry in the array determines whether a given neighbor is used for connectivity analyses. For example, connectivity = trues(3,3) would use 8-connectivity and test all pixels that touch the current one, even the corners.region = [1,3] would not test neighbors along dimension 2 for connectivity. This corresponds to just the nearest neighbors, i.e., 4-connectivity in 2d and 6-connectivity in 3d.The default is region = 1:ndims(A).
The output label is an integer array, where 0 is used for background pixels, and each connected region gets a different integer index.
colorset = distinguishable_colors(maximum(labeled_components)+1, colorant"black", dropseed = false)
labeled_components
480×854 Matrix{Int64}:
0 0 0 0 0 0 0 0 0 0 0 0 0 … 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 … 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 … 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249 249 249 249
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 249 249 249
⋮ ⋮ ⋮ ⋱ ⋮
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
coloredmask = map( n -> colorset[n+1], labeled_components)
maximum(labeled_components)
261
components_location = ImageTracking.locate_components(labeled_components)
min_area = 1
filter!(p -> p.area ≥ min_area, components_location)
length(components_location)
261
components_location = ImageTracking.locate_components(labeled_components)
min_area = 10
filter!(p -> p.area ≥ min_area, components_location)
length(components_location)
18
vd54 = copy(vd[54])
for p in components_location
draw!(vd54, Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
end
vd54
components_location = ImageTracking.locate_components(labeled_components)
min_area = 100
filter!(p -> p.area ≥ min_area, components_location)
3-element Vector{Main.ImageTracking.Blob}:
Blob Object
Centroid coordinates: (x, y) = (670.1450079239303, 387.1695721077655)
Span values: (xspan, yspan) = (640:701, 360:413)
Occupied: area = 2524
Blob Object
Centroid coordinates: (x, y) = (842.4875, 20.16875)
Span values: (xspan, yspan) = (837:849, 1:48)
Occupied: area = 160
Blob Object
Centroid coordinates: (x, y) = (850.0762711864406, 7.669491525423729)
Span values: (xspan, yspan) = (844:854, 1:17)
Occupied: area = 118
vd54 = copy(vd[54])
for p in components_location
draw!(vd54, Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
end
vd54
components_location = ImageTracking.locate_components(labeled_components)
min_area = 200
filter!(p -> p.area ≥ min_area, components_location)
1-element Vector{Main.ImageTracking.Blob}:
Blob Object
Centroid coordinates: (x, y) = (670.1450079239303, 387.1695721077655)
Span values: (xspan, yspan) = (640:701, 360:413)
Occupied: area = 2524
p = components_location[1]
draw(vd[54], Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
filename = "img/pendulo_70cm_1_reduzido_curto.mp4"
tracks_70cm1 = ImageTracking.find_tracks(filename)
Processing 150 frames ...100%|██████████████████████████| Time: 0:00:47
1-element Vector{Main.ImageTracking.Track}:
Tracked blob with framespan nspan = 13:150
ImageTracking.tracks_on_video(filename, tracks_70cm1)
┌ Info: Video saved with 150 frames and 1 tracks └ @ Main.ImageTracking /Users/rrosa/Documents/git_repositories/modelagem_matematica/notas_de_aula/scripts/ImageTracking.jl:150
DisplayMovie.display_movie("img/pendulo_70cm_1_reduzido_curto_tracked.mp4")
tracks_70cm1
1-element Vector{Main.ImageTracking.Track}:
Tracked blob with framespan nspan = 13:150
tracks_70cm1[1].nspan
13:150
tracks_70cm1[1].path
138-element Vector{Main.ImageTracking.Blob}:
Blob Object
Centroid coordinates: (x, y) = (848.7265500794913, 209.04451510333863)
Span values: (xspan, yspan) = (839:854, 175:241)
Occupied: area = 629
Blob Object
Centroid coordinates: (x, y) = (847.3719862227325, 212.89322617680827)
Span values: (xspan, yspan) = (836:854, 176:248)
Occupied: area = 871
Blob Object
Centroid coordinates: (x, y) = (846.1327913279133, 216.11020776874435)
Span values: (xspan, yspan) = (833:854, 177:254)
Occupied: area = 1107
Blob Object
Centroid coordinates: (x, y) = (844.8869047619048, 218.23511904761904)
Span values: (xspan, yspan) = (830:854, 177:256)
Occupied: area = 1344
Blob Object
Centroid coordinates: (x, y) = (843.5773130544993, 220.33523447401774)
Span values: (xspan, yspan) = (828:854, 177:258)
Occupied: area = 1578
Blob Object
Centroid coordinates: (x, y) = (842.7289227166276, 221.2980093676815)
Span values: (xspan, yspan) = (826:854, 177:258)
Occupied: area = 1708
Blob Object
Centroid coordinates: (x, y) = (842.6964285714286, 220.60059523809525)
Span values: (xspan, yspan) = (826:854, 178:258)
Occupied: area = 1680
Blob Object
Centroid coordinates: (x, y) = (842.8163514338011, 220.41122635753507)
Span values: (xspan, yspan) = (826:854, 178:258)
Occupied: area = 1639
Blob Object
Centroid coordinates: (x, y) = (842.8624923265808, 220.86126457949663)
Span values: (xspan, yspan) = (826:854, 179:259)
Occupied: area = 1629
Blob Object
Centroid coordinates: (x, y) = (842.980625, 221.265625)
Span values: (xspan, yspan) = (826:854, 179:258)
Occupied: area = 1600
Blob Object
Centroid coordinates: (x, y) = (842.8601615910503, 221.1224362958359)
Span values: (xspan, yspan) = (826:854, 179:258)
Occupied: area = 1609
Blob Object
Centroid coordinates: (x, y) = (842.5529622980251, 221.25852782764812)
Span values: (xspan, yspan) = (825:854, 179:259)
Occupied: area = 1671
Blob Object
Centroid coordinates: (x, y) = (842.5005966587112, 221.49224343675417)
Span values: (xspan, yspan) = (825:854, 179:259)
Occupied: area = 1676
â‹®
Blob Object
Centroid coordinates: (x, y) = (537.5756442227764, 423.214463840399)
Span values: (xspan, yspan) = (507:569, 399:446)
Occupied: area = 2406
Blob Object
Centroid coordinates: (x, y) = (586.1208315288003, 414.7986141186661)
Span values: (xspan, yspan) = (556:617, 391:438)
Occupied: area = 2309
Blob Object
Centroid coordinates: (x, y) = (631.0521201413428, 402.50132508833923)
Span values: (xspan, yspan) = (602:661, 379:426)
Occupied: area = 2264
Blob Object
Centroid coordinates: (x, y) = (671.4093648585999, 387.42280945757994)
Span values: (xspan, yspan) = (644:701, 364:411)
Occupied: area = 2157
Blob Object
Centroid coordinates: (x, y) = (707.1664247701983, 370.4446057087566)
Span values: (xspan, yspan) = (681:735, 348:394)
Occupied: area = 2067
Blob Object
Centroid coordinates: (x, y) = (737.0905852417303, 352.37150127226465)
Span values: (xspan, yspan) = (712:764, 330:376)
Occupied: area = 1965
Blob Object
Centroid coordinates: (x, y) = (762.6097691894794, 334.12238325281805)
Span values: (xspan, yspan) = (738:788, 312:358)
Occupied: area = 1863
Blob Object
Centroid coordinates: (x, y) = (783.1237172177879, 317.3369441277081)
Span values: (xspan, yspan) = (759:808, 296:340)
Occupied: area = 1754
Blob Object
Centroid coordinates: (x, y) = (798.8568935427575, 302.46364165212333)
Span values: (xspan, yspan) = (775:823, 281:325)
Occupied: area = 1719
Blob Object
Centroid coordinates: (x, y) = (810.6108343711084, 290.052303860523)
Span values: (xspan, yspan) = (787:834, 269:311)
Occupied: area = 1606
Blob Object
Centroid coordinates: (x, y) = (818.3419857235561, 280.7988319273199)
Span values: (xspan, yspan) = (796:841, 260:302)
Occupied: area = 1541
Blob Object
Centroid coordinates: (x, y) = (822.5608930987821, 274.98240866035184)
Span values: (xspan, yspan) = (800:845, 255:295)
Occupied: area = 1478
tracks_70cm1[1].path[1]
Blob Object Centroid coordinates: (x, y) = (848.7265500794913, 209.04451510333863) Span values: (xspan, yspan) = (839:854, 175:241) Occupied: area = 629
plt = plot(title = "Paths of $(length(tracks_70cm1)) track(s)", titlefont = 10,
xlims = (1,854), ylims = (1, 480),
size = (854, 480), aspect = :equal, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, getfield.(tr.path, :x), 481 .- getfield.(tr.path, :y), label="track $n")
end
display(plt)
plt = plot(title = "Coordinate x of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topleft)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, tr.nspan, getfield.(tr.path, :x), label="track $n")
end
display(plt)
plt = plot(title = "Coordinate y of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, tr.nspan, - getfield.(tr.path, :y), label="track $n")
end
display(plt)
size(vd[1])
(480, 854)
data_x = getfield.(tracks_70cm1[1].path, :x)
data_y = size(vd[1], 1) + 1 .- getfield.(tracks_70cm1[1].path, :y)
nothing
n_minima = Int[]
for j in 40:length(data_x)-1
if data_x[j-1] > data_x[j] < data_x[j+1]
push!(n_minima, j + first(tracks_70cm1[1].nspan) - 1)
end
end
n_minima
2-element Vector{Int64}:
72
125
plt = plot(title = "Coordinate x of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topleft)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, tr.nspan, getfield.(tr.path, :x), label="track $n")
end
vline!(plt, n_minima, label = "minima")
display(plt)
(n_minima[2:end] - n_minima[1:end-1]) / 29.5
1-element Vector{Float64}:
1.7966101694915255
g = 9.8 # m/s
l = 0.7 # m
T_p = 2π * √(l / g)
1.679251908362714
function medias(a::Real, b::Real)
(a ≥ 0 && b ≥ 0) || throw(ArgumentError("arguments must be nonnegative"))
ma = (a + b)/2
mg = sqrt(a * b)
return ma, mg
end
function agm(a::Real, b::Real; tol::Real = 1e-10, maxiter::Int = 100)
tol > 0 || throw(ArgumentError("tolerance must be positive"))
maxiter > 0 || throw(ArgumentError("maximum number of iterations must be positive"))
y1, y2 = a, b
n = 0
while abs(y1 - y2) > tol && n < maxiter
y1, y2 = medias(y1, y2)
n += 1
end
return (y1 + y2)/2, abs(y1 - y2), n
end
agm (generic function with 1 method)
θ₀ = π / 4
T_p / agm(1, cos(θ₀ / 2))[1]
1.7463772212095845
error(u) = sum(
abs2,
(data_x[50:end] .- u[1]).^2 + (data_y[50:end] .- u[2]).^2 .- u[3]^2
) / length(data_x)
error (generic function with 1 method)
center_point = Optim.optimize(error, [400.0, 800.0, 500])
* Status: success
* Candidate solution
Final objective value: 2.179085e+06
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 152
f(x) calls: 280
Optim.minimum(center_point)
2.1790845446313173e6
center_nx, center_ny, radius_n = Optim.minimizer(center_point)
3-element Vector{Float64}:
466.34634496986183
558.3533645572189
505.51516552719033
plt = plot(title = "Paths of $(length(tracks_70cm1)) track(s)", titlefont = 10,
xlims = (1,854), ylims = (1, 600),
size = (854, 600), aspect = :equal, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, data_x, data_y, label="track $n")
end
scatter!(plt, (center_nx, center_ny), markersize = 10)
plot!(plt, 1:854, nx -> center_ny - √(radius_n^2 - (nx - center_nx)^2), linestyle = :dash)
display(plt)
scale = radius_n / l # points per meter
722.1645221817005
data_scaled_x = (data_x .- center_nx) / scale
data_scaled_y = (data_y .- center_ny) / scale
138-element Vector{Float64}:
-0.3965825942201275
-0.4019120045625831
-0.4063666426583695
-0.40930906258292526
-0.41221714704552687
-0.41355032648607204
-0.412584598998547
-0.4123223749834595
-0.4129455546165464
-0.41350548300970785
-0.41330720588619096
-0.4134956553705288
-0.413819287454254
â‹®
-0.6931492935783854
-0.6814956475417342
-0.6644672715240699
-0.6435876586829319
-0.6200774982868892
-0.5950512004262687
-0.569781171978575
-0.5465379377714996
-0.5259424889246751
-0.5087561866204509
-0.49594266331796594
-0.48788851071379713
function dudt_pendulum!(du, u, p, t)
θ, ω = u
l, g, α = p
du[1] = ω
du[2] = - (g / l) * sin(θ) - α * ω
end
dudt_pendulum! (generic function with 1 method)
g = 9.8 # 9,8 m/s^2
l = 0.7 # 70 cm
α = 0
p = [l, g, α]
θ₀ = π/6
ω₀ = 0.0
u₀ = [θ₀, ω₀]
tspan = (0.0, 5.0) # até 5 segundos
prob_pendulum = ODEProblem(dudt_pendulum!, uâ‚€, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true timespan: (0.0, 5.0) u0: 2-element Vector{Float64}: 0.5235987755982988 0.0
sol = solve(prob_pendulum, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 32-element Vector{Float64}:
0.0
0.00014258482589892232
0.0015684330848881453
0.015826915674780374
0.05520287938356151
0.11869948603402829
0.2002131852019463
0.305914875953137
0.42890164394056673
0.5569187047083594
0.7013506994917333
0.872697562411765
1.0408497431878225
â‹®
2.6895148788362473
2.8890581891668203
3.1312540989213207
3.3360859425870872
3.579811053395161
3.7899434234522027
4.033677808420313
4.248777258478439
4.489724412426027
4.7129429751996135
4.9487942885038185
5.0
u: 32-element Vector{Vector{Float64}}:
[0.5235987755982988, 0.0]
[0.5235987044417862, -0.000998093740288535]
[0.5235901656815024, -0.01097897701799843]
[0.52272227807514, -0.1107323314966049]
[0.512965861754178, -0.3840406673395072]
[0.47498717025387077, -0.8072478892241957]
[0.388974682500463, -1.2882714556007608]
[0.22681168930963705, -1.74188417532362]
[-0.0034489375342210073, -1.936781390502851]
[-0.2416589452693402, -1.7139641695093224]
[-0.44371540806352233, -1.0197244000033383]
[-0.5224044417637903, 0.12912618721927982]
[-0.40599595484510004, 1.214547948859705]
â‹®
[-0.468181398562196, 0.8593957023872523]
[-0.19043787415408014, 1.8015157078868282]
[0.26158347131616566, 1.6731851139158112]
[0.5009792313597429, 0.5577249920462498]
[0.4331262533935047, -1.0800627594148071]
[0.10380928939438962, -1.8977491548843235]
[-0.337721682376994, -1.4733201977988184]
[-0.5219373906746808, -0.15523522675671586]
[-0.3642262071170868, 1.384044759406596]
[0.028820747076412908, 1.9341611787299822]
[0.4182220223013869, 1.1572024377115584]
[0.4697150849008621, 0.8484644252996695]
display(plot(sol, label = ["θ(t)" "ω(t)"]))
display(plot(sol, vars = 1, label = "θ(t)"))
n0 = 34
plt = plot(sol, vars = 1, label = "θ(t)")
scatter!(plt, (0:length(data_x)-n0) ./ 29.5, data_scaled_x[n0:end], label="tracked")
display(plt)
df = CSV.read("data/pendulo_70cm_reduzido.csv", DataFrame)
481 rows × 3 columns
| nf | nx | ny | |
|---|---|---|---|
| Int64 | Float64 | Float64 | |
| 1 | 13 | 848.727 | 209.073 |
| 2 | 14 | 847.362 | 212.843 |
| 3 | 15 | 846.127 | 216.006 |
| 4 | 16 | 844.89 | 218.281 |
| 5 | 17 | 843.577 | 220.369 |
| 6 | 18 | 842.701 | 221.178 |
| 7 | 19 | 842.682 | 220.634 |
| 8 | 20 | 842.833 | 220.468 |
| 9 | 21 | 842.866 | 220.873 |
| 10 | 22 | 843.008 | 221.232 |
| 11 | 23 | 842.875 | 221.118 |
| 12 | 24 | 842.552 | 221.235 |
| 13 | 25 | 842.523 | 221.488 |
| 14 | 26 | 842.712 | 221.488 |
| 15 | 27 | 842.606 | 221.49 |
| 16 | 28 | 842.406 | 221.461 |
| 17 | 29 | 842.366 | 221.539 |
| 18 | 30 | 842.363 | 221.584 |
| 19 | 31 | 842.25 | 221.574 |
| 20 | 32 | 842.034 | 221.523 |
| 21 | 33 | 841.969 | 221.37 |
| 22 | 34 | 841.981 | 221.283 |
| 23 | 35 | 841.954 | 221.323 |
| 24 | 36 | 842.037 | 221.198 |
| 25 | 37 | 842.135 | 221.06 |
| 26 | 38 | 842.189 | 221.251 |
| 27 | 39 | 842.266 | 221.351 |
| 28 | 40 | 842.162 | 221.446 |
| 29 | 41 | 841.951 | 221.716 |
| 30 | 42 | 841.715 | 221.619 |
| ⋮ | ⋮ | ⋮ | ⋮ |
plt = plot(title = "Coordinate x", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, df.nx, label = nothing)
display(plt)
plt = plot(title = "Coordinate y", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, size(vd[1], 1) .+ 1 .- df.ny, label = nothing)
display(plt)
data2_scaled_t = (df.nf[n0:end] .- df.nf[n0]) / 29.97
data2_scaled_x = (df.nx[n0:end] .- center_nx) / scale
data2_scaled_y = (size(vd[1], 1) .+ 1 .- df.ny[n0:end] .- center_ny ) / scale
display(scatter(data2_scaled_t, data2_scaled_x))
display(scatter(data2_scaled_t, data2_scaled_y))
data_x2 = getfield.(tracks2[1].path, :x)
n_minima2 = Int[]
for j in 40:length(data_x2)-1
if data_x2[j-1] > data_x2[j] < data_x2[j+1]
push!(n_minima2, j + first(tracks2[1].nspan) - 1)
end
end
n_minima2
n_minima2 = Int[]
for j in 40:length(df.nx)-1
if df.nx[j-1] > df.nx[j] < df.nx[j+1]
push!(n_minima2, j + df.nf[1] - 1)
end
end
n_minima2
8-element Vector{Int64}:
72
125
178
231
283
336
389
441
(n_minima2[2:end] - n_minima2[1:end-1])/30
7-element Vector{Float64}:
1.7666666666666666
1.7666666666666666
1.7666666666666666
1.7333333333333334
1.7666666666666666
1.7666666666666666
1.7333333333333334
sum(n_minima2[2:end] - n_minima2[1:end-1])/7/30
1.7571428571428571
plt = plot(title = "Coordinate x", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, df.nx, label = "data")
vline!(plt, n_minima2, label = "minima")
display(plt)
g = 9.8 # 9,8 m/s^2
l = 0.7 # 70 cm
α = 0
p = [l, g, α]
θ₀ = π/6
ω₀ = 0.0
u₀ = [θ₀, ω₀]
tspan2 = (0.0, 16.45)
prob_pendulum2 = ODEProblem(dudt_pendulum!, uâ‚€, tspan2, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true timespan: (0.0, 16.45) u0: 2-element Vector{Float64}: 0.5235987755982988 0.0
VideoIO.get_duration("img/pendulo_70cm_1_reduzido.mov")
16.45
θ₀ = π / 6
α = 0.0 #-0.1
remake(prob_pendulum2; u0 = [θ₀, 0.0], p = [l, g, α])
sol2 = solve(prob_pendulum2, Tsit5(); saveat = 1/29.87)
data2_scaled_x[1:10] .- sol2(data2_scaled_t)[1,1:10]
#length(sol[1,:])
#sol(data2_scaled_t)[1,:]
sol(data2_scaled_t)
t: 448-element Vector{Float64}:
0.0
0.033366700033366704
0.06673340006673341
0.1001001001001001
0.13346680013346682
0.16683350016683351
0.2002002002002002
0.2335669002335669
0.26693360026693363
0.3003003003003003
0.33366700033366703
0.3670337003670337
0.4004004004004004
â‹®
14.547881214547882
14.58124791458125
14.614614614614615
14.647981314647982
14.68134801468135
14.714714714714715
14.748081414748082
14.78144811478145
14.814814814814815
14.848181514848182
14.88154821488155
14.914914914914915
u: 448-element Vector{Vector{Float64}}:
[0.5235987755982988, 0.0]
[0.5197064819391034, -0.2330414413613375]
[0.5080822001687944, -0.4629302691579599]
[0.488883698677749, -0.6865155884298596]
[0.46237373952364935, -0.9006563907398047]
[0.4289196227578753, -1.102237149997204]
[0.38899141026001377, -1.288202512164673]
[0.3431598607271923, -1.4556077279142619]
[0.29209180679967006, -1.6016605695396418]
[0.23654135225424, -1.7238352335082683]
[0.17734352255250205, -1.819932863153441]
[0.11540186183275535, -1.88814487747045]
[0.05166798968064211, -1.9271578012356183]
â‹®
[26920.697182912332, 33556.85056390745]
[27306.360007144147, 33989.425031612314]
[27696.14382775103, 34426.118397211016]
[28090.07784082385, 34866.95635621334]
[28488.191345596126, 35311.96468510579]
[28890.51374518916, 35761.16921555201]
[29297.07454454987, 36214.59588839963]
[29707.903351789308, 36672.27070221071]
[30123.029879305846, 37134.21973321311]
[30542.483937377587, 37600.4691636604]
[30966.29544585134, 38071.04522386799]
[31394.494424368553, 38545.97422357512]
plot(data2_scaled_t, sol(data2_scaled_t)[1,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
function ddu_pendulum!(ddu, du, u, p, t)
θ = u[1]
ω = du[1]
l, g, α = p
ddu[1] = - (g / l) * sin(θ) - α * ω
end
ddu_pendulum! (generic function with 1 method)
θ₀ = π / 6
ω₀ = 0.0
α = 0.05
l = 0.74
p = [l, g, α]
prob_pendulum_2nd = SecondOrderODEProblem(ddu_pendulum!, [ω₀], [θ₀], tspan2, p)
sol_2nd = solve(prob_pendulum_2nd)
sol_2nd.retcode
:Success
plot(data2_scaled_t, sol_2nd(data2_scaled_t)[2,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
prob_rmk = remake(prob_pendulum_2nd; u0 = [θ₀], p = [l, g, α])
sol = solve(prob_rmk)
type Array has no field x
Stacktrace:
[1] getproperty(x::Vector{Float64}, f::Symbol)
@ Base ./Base.jl:33
[2] (::DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing})(du::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}, t::Float64)
@ SciMLBase ~/.julia/packages/SciMLBase/UIp7W/src/scimlfunctions.jl:354
[3] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}}, Tsit5, OrdinaryDiffEq.InterpolationData{DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/perform_step/low_order_rk_perform_step.jl:623
[4] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:default_set, :second_time), Tuple{Bool, Bool}}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:455
[5] #__solve#404
@ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:4 [inlined]
[6] __solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}}, ::Nothing; default_set::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:second_time,), Tuple{Bool}}})
@ DifferentialEquations ~/.julia/packages/DifferentialEquations/7WPQg/src/default_solve.jl:8
[7] #__solve#57
@ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:291 [inlined]
[8] __solve
@ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:278 [inlined]
[9] #solve_call#42
@ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:61 [inlined]
[10] solve_call
@ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:48 [inlined]
[11] #solve_up#44
@ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:90 [inlined]
[12] solve_up
@ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:78 [inlined]
[13] #solve#43
@ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:73 [inlined]
[14] solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:68
[15] top-level scope
@ In[137]:2
[16] eval
@ ./boot.jl:360 [inlined]
[17] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1116
function objective(β)
g = 9.8
ω₀ = 0.0
θ₀, l, α = β
p = [l, g, α]
tspan2 = (0.0, 16.45)
prob = remake(prob_pendulum_2nd, u0 = [θ₀, 0.0], p = [l, g, α])
#sol = solve(prob)
prob = SecondOrderODEProblem(ddu_pendulum!, [ω₀], [θ₀], tspan2, p)
sol = solve(prob)
return mse(data2_scaled_x, sol(data2_scaled_t)[2,:])
end
objective (generic function with 1 method)
res = Optim.optimize(objective, [Ï€ / 3, 0.7, 0.0])
* Status: success
* Candidate solution
Final objective value: 1.313365e-03
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 69
f(x) calls: 137
θ₀, l, α = Optim.minimizer(res)
3-element Vector{Float64}:
0.5557462725091609
0.7421821761291346
0.03752799438015029
π / θ₀
5.652926180513441
θ₀ - π / 5.6529
-2.57385107782504e-6
remake(prob_pendulum_2nd, u0 = [θ₀, 0.0], p = [l, g, α])
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true timespan: (0.0, 16.45) u0: 2-element Vector{Float64}: 0.5557462725091609 0.0
sol_2nd = solve(prob_pendulum_2nd)
sol_2nd.retcode
:Success
plot(data2_scaled_t, sol_2nd(data2_scaled_t)[2,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")